// B2Function.cpp: implementation of the B2Function class.
//
//////////////////////////////////////////////////////////////////////
#include <iomanip>
#include <fstream>
#include <string>
#include "B2_function.h"
#include "DataHolder.h"


//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

B2Function::B2Function(int i_NInput)
{
  NInput = i_NInput;
  nCall = 0;
  nMult = 1;

  // Create 3 arrays
  newA(mCosArray, NInput);          // stores the cosine values
  newA(mSinArray, NInput);          // stores the sine values
  newA(mDerivative, NInput);        // a dummy derivative used in constructing the derivative

  Result.Size(NInput);
}


B2Function::~B2Function()
{
  deleteA(mCosArray);
  deleteA(mSinArray);
  deleteA(mDerivative); 
}

// takes a point and if it is beyong boundary move the point
// into the boundary box
void B2Function::checkBoundary(vector<double> &Input)
{
  for (unsigned int i = 0; i < Input.size(); i++)
  {
    // reflective boundary
    /*while( fabs(Input[i]) >= 100.0 )
    {
      if ( Input[i] >= 100.0 )
      {
        Input[i] = 100.0 - fabs(Input[i]-100.0);
      }
      else if ( Input[i] <= -100.0 )
      {
        Input[i] = -100.0 + fabs(-100.0-Input[i]);
      }
    }*/

    // periodic boundary
    if ( Input[i] >= 1.0 )
    {
      Input[i] = -1.0 + fmod(Input[i],1.0);
    }
    else if ( Input[i] <= -1.0 )
    {
      Input[i] = 1.0 + fmod(Input[i],1.0);
    }
  }
}

// takes a point and if it is beyong boundary move the point
// into the boundary box
void B2Function::checkBoundary(double *Input)
{
  for (int i = 1; i <= NInput; i++)
  {
    // periodic boundary
    if ( Input[i] >= 1.0 )
    {
      Input[i] = -1.0 + fmod(Input[i],1.0);
    }
    else if ( Input[i] <= -1.0 )
    {
      Input[i] = 1.0 + fmod(Input[i],1.0);
    }
  }
}

//////////////////////////////////////////////////////////////////////
// functions
//////////////////////////////////////////////////////////////////////

double B2Function::Value(double *Input)
{
  // first check that point is not beyong any boundaries
  checkBoundary(Input);

  // Initialise vars answer and danswer (both local)
  double answer = 0.0;
  double danswer = 0.0;
  double arg1 = 3 * PI;
  double arg2 = 4 * PI;

  answer += Input[1]*Input[1];
  danswer += Input[2]*Input[2];
  danswer *= 2;
  answer += danswer;
  answer -= 0.3*cos(arg1*Input[1]);
  answer -= 0.4*cos(arg2*Input[2]);
  answer += 0.7;

  Result.Input(answer, Input);

  nCall++;

  return (answer); 
}

// calculates/builds the function value
double B2Function::Value(vector<double> &Input)
{
  // first check that point is not beyong any boundaries
  checkBoundary(Input);

  // Initialise vars answer and danswer (both local)
  double answer = 0.0;
  double danswer = 0.0;
  double arg1 = 3 * PI;
  double arg2 = 4 * PI;

  answer += Input[1]*Input[1];
  danswer += Input[2]*Input[2];
  danswer *= 2;
  answer += danswer;
  answer -= 0.3*cos(arg1*Input[1]);
  answer -= 0.4*cos(arg2*Input[2]);
  answer += 0.7;

  Result.Input(answer, Input);

  nCall++;

  return (answer); 
}

// calculates/builds the derivative (used for CGM and HMC)
void B2Function::BuildDerivative(double *Input, double *Derivative)
{
  // NInput = dimension of function
  double mDer1 = 0.9 * PI;
  double mDer2 = 1.6 * PI;
  double arg1 = 3 * PI;
  double arg2 = 4 * PI;

  // Constructing equation
  for(int i=1; i<=NInput; i++)
  {
    Derivative[i] = 0.0;
    mDerivative[i] = 0.0;
  }

  Derivative[1] += 2 * Input[1];
  Derivative[1] += mDer1 * sin(arg1 * Input[1]);
  Derivative[2] += 4 * Input[2];
  Derivative[2] += mDer2 * sin(arg2 * Input[2]);

  for(int j = 1; j<=NInput; j++)
  {
    mDerivative[j] = Derivative[j];
  }
}

// used to test if the function is working properly
void B2Function::test_function()
{
  ofstream out("output/test_of_b2.txt");
  
  out << "NUMBER_OF_DIMENSIONS: " << NInput << endl 
    << "      Subscript[x,k]      |     Value    |         Derivative[k]" << endl << endl;

  int NRun = 3;
  double *mInput;
  double *mDeriv;
  newA(mInput, NInput);
  newA(mDeriv, NInput);

  for(int run=1; run<=NRun; run++)
  {
    for(int n=1; n<=NInput; n++)
    {
      mInput[n]=(2.0*ran1())-1.0;
    }

    for(int k=1; k<=NInput; k++)
    {
      out << setw(11) << mInput[k];
    }
        
    out << setw(16) << Value(mInput);

    BuildDerivative(mInput, mDeriv);

    for(int l=1; l<=NInput; l++)
    {
      out << setw(16) << mDerivative[l];
    }

    out << "\n\n";
    
  }

  deleteA(mInput);
  deleteA(mDeriv);
  out.close();
}

// used to print the solution (this is the old version, the new version is below)
void B2Function::print_sol(double *solution, double funcVal, char filename[])
{
  ofstream out(filename);	

  out << "Number of dimensions: " << NInput << endl;
  out << "Final B2 value = " << funcVal << endl;
  out << "Solution coordinates:" << endl;

  for (int i = 1; i <= NInput; i++)
  {
    out << solution[i] << "  ";
  }
 
  out << endl;

  out.close();
}

// used to print the solution to a file
void B2Function::print_sol(char filename[])
{
  ofstream out(filename);	

  // getBestValue() returns the best value stored in the tracker
  out << "Number of dimensions: " << NInput << endl;
  out << "Final B2 value = " << Result.getBestValue() << endl;
  out << "Solution coordinates:" << endl;

  // prints out the solution coordinates to a file
  Result.printSolution(out);

  out.close();
}


// distance from know global minimum which is located at the origin
double B2Function::distanceFromGlobalMinimum(double *solution)
{
  double toReturn = 0.0;

  for (int i = 1; i <= NInput; i++)
  {
    toReturn += solution[i]*solution[i];
  }

  double toReturnS = sqrt(toReturn);
  return toReturnS;
}